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f*"*) , A new iterative method is developed to numerically calculate the periodic, matched beam envelope 

C^l ■ solution of the coupled Kapchinskij-Vladimirskij (KV) equations describing the transverse evolution 

of a beam in a periodic, linear focusing lattice of arbitrary complexity. Implementation of the method 
^ , is straightforward. It is highly convergent and can be applied to all usual parameterizations of the 

matched envelope solutions. The method is applicable to all classes of linear focusing lattices without 
skew couplings, and also applies to all physically achievable system parameters - including where 
the matched beam envelope is strongly unstable. Example applications are presented for periodic 
solenoidal and quadrupole focusing lattices. Convergence properties are summarized over a wide 
range of system parameters. 
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C0 ■ I. INTRODUCTION 

o ■ 

The Kapchinskij-Vladimirskij (KV) envelope equations Q, 0, are often employed as a simple description of the 
f-H , transverse evolution of intense ion beams. The equations are coupled ordinary differential equations that describe the 
O ■ evolution of the beam edge (or rms radii) in response to applied linear focusing forces of the lattice and defocusing 
— i \ forces resulting from beam space-charge and transverse phase-space area (emittances). Although the KV envelope 

| ■ equations are only fully Vlasov consistent with the singular KV distribution, the equations can be applied to describe 

the low-order evolution of a real distribution of beam particles when the variation of the statistical beam emittances is 
negligible or sufficiently slowjsj- Large nonlinear fields that can be produced by non-ideal applied focusing elements, 
(/*■) ' nonuniform beam space-charge, and species contamination (electron cloud effects, etc.) drive deviations from the KV 
i-H model. Such effects arc suppressed to the extent possible in most practical designs, rendering the KV model widely 
CN ■ applicable. 

The matched solution of the KV envelope equations is the solution with the same periodicity as the focusing 
latticeElllQ. The matched beam envelope is important because it is believed to be the most radially compact solution 
supported by a periodic linear focusing channel^- Matched envelopes are typically calculated as a first step in the 
design of practical transport lattices and for use in initializing more detailed beam simulations to evaluate machine 
<fh ' performance 2]. The matched envelope solution is typically calculated by numerically integrating trial solutions of 
the KV equations from assumed initial conditions over one lattice period and searching for the four initial envelope 
^h* 1 coordinates and angles that generate the solution with the periodicity of the lattice^, 0, 0- An optimal formulation 
of the conventional root finding procedure for envelope matching has been presented by RyneQ. Conventional root 
• • ■ finding procedures for matching can be surprisingly problematic even for relatively simple focusing lattices. Variations 
in initial conditions can lead to many inflection points in the envelope functions at the end of the lattice period. Thus 
initial guesses close to the actual values corresponding to the periodic solution are often necessary to employ standard 
root finding techniques. This is especially true for complicated focusing lattices with low degrees of symmetry and 
where focusing strengths (or equivalently, undepressed single particle phase advances) are large. For large focusing 
strength and strong space-charge intensity, the matched envelope solution can be unstable over a wide range of system 
parametersj^OI- Such instabilities can restrict the basin of attraction when standard numerical root finding methods 
are used to calculate the needed matching conditions — especially for certain classes of solution parameterizations. 

In this article we present a new iterative procedure to numerically calculate matched envelope solutions of the KV 
equations. The basis of this procedure is that the KV distribution of particle orbits internal to the beam must have 
a locus of turning points consistent with the beam edge (envelope). In the absence of beam space-charge, betatron 
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amplitudes calculated from the sine- and cosine-like principal orbits describing particles moving in the applied focusing 
fields of the lattice directly specify the matched beam envelope^. For finite beam space-charge, the principal orbits 
describing the betatron amplitudes and matched beam envelope cannot be a priori calculated because the defocusing 
forces from beam space-charge uniformly distributed within the (undetermined) beam envelope are unknown. In 
the iterative matching (IM) method, the relation between the betatron amplitudes and the particle orbits is viewed 
as a consistency equation. Starting from a simple trial envelope solution that accounts for both space-charge and 
and applied focusing forces in a general manner, the consistency condition is used to iteratively correct the envelope 
functions until converged matched envelope solutions are obtained that are consistent with particle orbits internal to 
the beam. 

The IM method offers superior performance and reliability in constructing matched envelopes over conventional 
root finding because the IM iterations are structured to reflect the periodicity of the actual matched solution rather 
than searching for parameters that lead to periodicity. The IM method works for all physically achievable system 
parameters (even in cases of envelope instability) and is most naturally expressed and rapidly convergent when relative 
beam space-charge strength is expressed in terms of the depressed particle phase advance. All other parameterizations 
of solutions (specified perveances and emittances, etc.) can also be carried out by simple extensions of the IM method 
rendering the approach fully general. The natural depressed phase advance parameterization is also useful when 
carrying out parametric studies because phase advances are the most relevant parameters for analysis of resonance- 
like effects central to charged particle dynamics in accelerators. The IM method provides a complement to recent 
analytical perturbation theories developed to construct matched beam envelopes in lattices with certain classes of 
symmetries |(J 0, 0, 0. I n contrast to these analytical theories, the IM method can be applied to arbitrary linear 
focusing lattices without skew couplings. The highly convergent iterative corrections of the IM method have the same 
form for all order iterations after seeding, rendering the method straightforward to code and apply to numerically 
generate accurate matched envelope solutions. 

The organization of this paper is the following. After a review of the KV envelope equations in Sec. [HJ various 
properties of matched envelope solutions and the continuous focusing limit arc analyzed in Sec. II I II These results 
are used in Sec. IIVI to formulate the IM method for calculation of matched solutions to the KV envelope equations. 
Example applications of the IM method are presented in Sec. to illustrate application and convergence properties 
of the method over a wide range of system parameters for a variety of systems. Concluding comments in Sec IVII 
summarize the advantages of the IM method over conventional techniques. 



II. THEORETICAL MODEL 



We consider an unbundled beam of particles of charge q and mass m coasting with axial relativistic factors /3b = const 
and 7b = 1/^/1 — 0^. In the KV model, the beam is propagating in a linear focusing lattice without skew couplings 
and has uniform charge density within an elliptical cross-section with principal radii r x and r y along the (transverse) 
x- and y-coordinate axes. When self-fields are included and image effects are neglected, the envelope radii consistent 
with the KV distribution evolve according to the so-called KV envelope equations 0,0, 

lO e 2 

r'>{s) + KjtfrjW - prj^M JTT = °- (1) 

J r x {s) + ry(s) rj(s) 

Here, primes denote derivatives with respect to the axial machine coordinate s, the subscript j ranges over both 
transverse coordinates x and y, the functions Kj(s) represent linear applied focusing forces of the transport lattice, 
Q = const is the dimcnsionlcss pcrvcance, and £j = const arc the rms edge emittances. Equations relating the 
functions kj to magnetic and/or electric fields of practical focusing elements are presented in Ref. 0]. The perveance 
provides a dimensionless measure of self-field defocusing forces internal to the beam0| and is defined as 



*2 o_, ™„3„,3/33 ' ( 2 ) 



Here, I is the constant beam current, c is the speed of light in vacuo, and en is the permittivity of free space. The 
perveance Q can be thought of as a scaled measure of space-charge strengthQ- The rms edge emittances £j provide 
a statistical measure of beam phase-space area projections in the x-x' and y-y' planes 

When the emittances are constant (e, = const), the KV envelope equations Q are consistent with the Vlasov 
equation only for the KV distributionjll Il0l| . which is a singular function of Courant-Snyder invariants. This singular 
structure can lead to unphysical instabilities within the Vlasov model[Tl|. However, the KV envelope equations can 
be applied to physical (smooth) distributions in an rms equivalent beam sense 0, with the envelope radii and the 
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omittances defined by statistical averages of the physical distribution as 

r, = 2v^), r y = 2v/(7), (3) 



and 



2w I2\ I /\21 1 / 2 



e x = 4 [{x-){x-j - {xx, 
e y = 4[(y 2 )(y' 2 )-(yyT] 1/2 - 

Here, (• • • ) denotes a transverse statistical average over the beam distribution function, and for notational simplicity, 
we have assumed zero centroid offset (e.g., (x) =0). In this rms equivalent sense, the omittances Ej will generally 
evolve in s. If this variation has negligible effect on the Tj, then the KV envelope equations can be applied with 
Ej = const to reliably model practical machines. This must generally be verified a posteriori with simulations of the 
full distribution. 

For appropriate choices of the lattice focusing functions Kj(s), Eq. Q can be employed to model a wide range 
of transport channels, including solcnoidal and quadrupole transport. For solenoidal transport, the equations must 
be interpreted in a rotating Larmor frame (see Appendix A of Ref. Q). In a periodic transport lattice, the Kj are 
periodic with fundamental lattice period L p , i.e., 

Kj(s + L p ) = Kj(s). (5) 

The beam envelope is said to be matched to the transport lattice when the envelope functions have the same periodicity 
as the lattice: 

r j (s + L p )=r j (s). (6) 

For specified focusing functions Kj(s), beam perveancc Q, and omittances Ej, the matching condition is equivalent to 
requiring that rj and r'j satisfy specific initial conditions at s = Sj when the envelope equations Q are integrated 
as an initial value problem. The required initial conditions generally vary with the phase of Sj in the lattice period 
(because the conditions vary with the local matched solution). In conventional procedures for envelope matching, 
needed initial conditions are typically found by numerical root finding starting from guessed seed valuesQ- This 
numerical matching can be especially problematic when: applied focusing strengths are large, the focusing lattice 
is complicated and devoid of symmetries that can reduce the dimensionality of the root finding, choices of solution 
parameters require extra constraints to effect, and where the matched beam envelope is unstable. 

The undepressed particle phase advance per lattice period <7oi_provides a dimensionless measure of the strength of 
the applied focusing functions Kj describing the periodic lattice^, ■ The o~Qj can be calculated from[j| 

coseroj = ^Tr M aj (s l + L p \si), (7) 

where M(y(s|sj) denotes the 2x2 single particle transfer matrix in the j-plane from axial coordinate s, to s. Explicitly, 
we have 

_ / C 0j (s\si) S j(s\si) 



where the C(y(s|sj) and iS"oj (s|sj) denote cosine-like and sine-like principal orbit functions satisfying 

FH j (s\s i ) + K j (s)F oj (s\s i ) = 0, (9) 

with F representing C or S with Coj subject to cosine-like initial (s = Si) conditions Coj(si\si) = 1 and CQj(si\si) = 0, 
and with Soj subject to sine- like initial conditions Soj(si\si) = and S' 0: j(si\si) = 1. Equation Q can be expressed 
in terms of Coj and Sqj as 

cosa-oj = -^[C j(si + L p \si) + S' 0j (si + L p \si)}. (10) 

The o~oj are independent of the particular value of Si used in the calculation of the principal functions. For some 
particular cases such as piecewise constant Kj the principal functions Foj can be calculated analytically. But, in 
general, the F j must be calculated numerically. In the absence of space-charge, the particle orbit is stable whenever 



particle o 

o-Qj < 180° and parametric bands of stability can also usually be found for a oj > 180° 0, 0, E3 • 

For a stable orbit, the 
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scale of the Kj (i.e., Kj — > ctKj with a = const setting the scale of the specified Kj) can always be regarded as being 
set by the cr j- I n this context, Eq. (|10|l is employed to fix the scale of the Kj in terms of agj and other parameters 
defining the Kj. Because there appears to be no advantage in using stronger focusing with aoj > 180° in terms of 
producing more radially compact matched envelopes [3I Il3j| . we will assume in all analysis that follows that the Kj are 
sufficiently weak to satisfy 00 < 180°. 

The formulation given above for calculation of the undepressed principal orbits Coj and Soj and the undepressed 
particle phase advances aoj can also be applied to calculate the depressed principal orbits Cj and Sj and the depressed 
phase advances Oj in the presence of uniform beam space-charge density for a particle moving within the matched 
KV beam envelopes. This is done by replacing 

Kj — > Kj - 7 r (11) 

in Eqs. @ and dropping the subscript 0s in Eqs. Q- 110(1 for notational clarity (i.e., Coj — > Cj and Soj — > Sj). 
Explicitly, the depressed principal functions satisfy 

w^i.ww- J;;*,, ^ (i2) 

with F representing C or S with Cj subject to Cj(si\si) = 1 and Cj(s,-|sj) = 0, and Sj subject to Sj(si\si) = and 
S'j(si\si) = 1, and the depressed phase advances satisfy 

cosa ] = \\ C ]i s i + L p\ s i) + S 'j( s i + L p\ s i)}- (13) 
For a stable orbit, it can be shown that the Oj can also be calculated from the matched envelope as UHJ 

ds 

(14) 



rj(s) 

This formula can also be applied to calculate aoj by using the matched envelope functions Tj calculated with Q = 0. 

Matched envelope solutions of Eqs. Q can be regarded as being determined by the focusing functions Kj, the 
perveance Q, and the emittances Sj. The lattice period L p is implicitly specified through the Kj. We will always 
regard the scale of the Kj as being set by the undepressed phase advances cr j through Eq. I|1U|) . For cr j < 180° 
there is no ambiguity in scale choice and the use of the eroj as parameters allows disparate classes of lattices to be 
analyzed in a common framework^. The depressed phase advances a x and a y can be employed to replace up to 
two of the three parameters Q, e x , and e y . Such replacements can be convenient, particularly when carrying out 
parametric surveys (for example, see Ref. 0) because crj/aoj is a dimensionless measure of space-charge strength 
satisfying < crj/croj < 1 with (Tj/aoj — * 1 representing a warm beam with negligible space-charge (i.e., Q — > 0, or 
Sj — * 00 for finite Q), and cfj/aoj — > representing a cold beam with maximum space-charge intensity (i.e., Ej — > 0). 
We will discuss calculation of matched beam envelopes for the useful parameterization cases listed in Table [I] In cases 
typical of linear accelerators the focusing functions have equal strength in the x- and y-planes giving <tq x = a^ y . In 
such plane symmetric cases we denote eroj = ctq. In practical situations where the focusing lattice and emittances 
arc both plane symmetric with eroj = Co an d Sj = e, then the depressed phase advance is also plane symmetric with 
Uj = a and parameterization cases 2 and 3 are identical. It is assumed that a unique matched envelope solution 
exists independent of the parameterization when the Kj are fully specified. There is no known proof of this conjecture, 
but numerical evidence suggests that it is correct for simple focusing lattices (i.e., simple Kj) when aoj < 180°. In 
typical experimental situations, note that transport lattices are fixed in geometry and excitations of focusing elements 
in the lattices can be individually adjusted. In the language adopted here, such lattices with different excitations in 
focusing elements (both overall scale and otherwise) correspond to different lattices described by different Kj with 
corresponding different matched envelopes. 

III. MATCHED ENVELOPE PROPERTIES 

In development of the IM method in Sec. IIVI we employ a consistency equation between depressed particle orbits 
within the beam and the matched envelope functions (|III Afl and use a continuous focusing description of the matched 
beam (IIII B|) to construct a seed iteration. Henceforth, we denote lattice period averages with overbars, i.e., for some 
quantity C( s )> 

C=/ ^C(s). (15) 
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TABLE I: Possible parameterizations of matched envelope solutions. 
Case Parameters 

Kj (cr j), Q, Ej 

1 Kj (ffOj), Q, CTj 

2 Kj (croj), Sj, and one of Oj 

3 Kj {&oj), Oj, and one of ej 



A. Consistency condition between particle orbits and the matched envelope 

We calculate nonlinear consistency conditions for the matched envelope functions Tj and the depressed principal 
orbit functions Cj and Sj as follows. First, the transfer matrix Mj of the depressed particle orbit in the j-plane is 
expressed in terms of betatron function- like formulation as[f| 



with, 



Here. 



M,-(s| Si ) 



C j( 8 \ s i) Sj(s\si) 



Cj(s\si) 
Sj(s\si) 

S'Js\ Si ) 



rj(si) 



cos Ai(j j(s) 



sin Aipj(s), 



rj(si)rj(s) . . 
-A±—LJ^-J- sm (s) , 



r'(s) r'J Si ) 



rj(si) rj(s) 



cos Ai/jj (s) 



rj{si)rj(s) 



sin Aipj(s), 



rj(si) 



cos Aijjj(s) 



r j( s i) r 'i( s ) . 



sin Atpj(s). 



(16) 



(17) 



A^(s) 



ds 



(18) 



is the change in betatron phase of the particle orbit from ,s = s, to s and the principal functions Cj and Sj are 
calculated including the linear space-charge term of the uniform density elliptical beam from Eq. 1)12(1 . Note that 
Tj = \Jej(3j can be used in Eqs. ((T7jl and ((TH|) to express the results more conventionally in terms of the betatron 
amplitude functions /3j describing linear orbits internal to the beam in the j-planeQ. These generalized betatron 
functions are periodic [i.e., (3j(s + L p ) = (3j(s) and include the transverse defocusing effects of uniformly distributed 
space-charge within the KV equilibrium envelope. Recognizing that Aifjj(si + L p ) = <jj [see Eq. ()14f) ] and that the 
matched envelope functions Tj have period L p gives 



r|(s) _ [M^js + L p \s) _ Sjis 



l p\ s ) 



sm a \ 



sm a j 



(19) 



Here, [Mj]i2 denotes the 1,2 component of the 2x2 matrix Mj and Uj can be equivalently calculated from cither 
Eq. (HSJ) or Eq. 1(11)1. 

Equation (|19f) can be applied to numerically calculate the consistency conditions for the matched envelope functions 
Tj on a discretized axial grid of s locations. As written, the principal orbit functions employed (i.e., the Cj and 
Sj) need to be independently calculated at each s-location on the grid through one lattice period. The fact that 
every period is the same can be applied to simplify the calculation. For any initial axial coordinate Si we have 
Mj(s + L p \s) = Mj(s + L p \si + L p ) ■ Mj(si + L p \s). Multiplying this equation from the right side by the identity 
matrix I = Mj(s|sj) • M.J 1 (s\si) where M" 1 is the inverse matrix and using Mj(s.; + L p \s) ■ M.j(s\si) = M.j(si + L p \si) 
gives 



Mj(s + L p \s) = Mj(s\ Si ) ■ Mj( Si + L p \ Si ) ■ M . 1 {s\s l ). 



(20) 
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Some straightforward algebra employing Eqs. i|16|) . (|19|> - and 1|20[1 . and the Wronskian (or symplectic) condition on 

mM 



Cjisisjs'Msi) - SMsiiCM*) = 1 ( 21 ) 



yields 



&00 



r ( S ) 5 ( S | Sl ) 



-j ^(si + Lp|a»)/ sincrj 
^-(sj + L p |si) 



r , , s cosgj - Cj(si + L p |sj) 



Equation (|22[1 explicitly shows that the linear principal functions Cj and Sj need only be calculated in s from some 
arbitrary initial point (sj) over one lattice period (to Sj + L p ) to calculate the consistency condition for the matched 
envelope functions rj(s), or equivalently, the betatron amplitude functions f3j(s) = r|(s)/£j. Equation 1221) can also 
be derived using Courant-Snydcr invariants of particle orbits within the beam. 

Equations and i122li form the foundation of an iterative numerical method developed in Sec. IIVI to calculate 
the matched beam envelope for any lattice. These equations express the intricate connection between the bundle of 
depressed particle orbits within the uniform density KV beam and the locus of maximum particle excursions defining 
the envelope functions rj. The method will be iterative because the consistent matched envelope functions rj arc 
necessary to integrate the linear differential equations for the depressed orbit principal functions Cj and Sj. However, 
in the limit Q —> 0, the principal functions do not depend on the rj and the matched envelope can be immediately 
calculated from the equations. Thus, the periodic zero-current matched beam envelope can be directly calculated 
using Eq. (|22|l in terms of the two independent, aperiodic linear orbits (i.e., Coj and Soj) integrated over one lattice 
period. 

Additional constraints on the matched envelope functions rj and/or betatron functions (3j are necessary to formulate 
the IM method for parameterizations where one or more of the parameters Q and ej need to be eliminated (see TablcQJ . 
Appropriate constraints can be derived by taking the period average of Eq. for a matched envelope, giving 



2Q^— -£■- = 0. (23) 

r X + Ty J Tj 



B. Continuous limit 



In the continuous focusing approximation, we take the lattice focusing functions Kj as constants set according to 

with the <7(y calculated from Eq. (|10|l consistent with the actual s- varying periodic focusing functions Kj. Then we 
replace rj — > TJ in the KV envelope equations Q and take rj = const to obtain the continuous limit envelope equation 



2*y^_ = j5 = _4=o. (25) 

L P J 3 r x +r y rj 3 

Equation l|25() provides an estimate of the lattice period average envelope radii rj in response to applied focusing forces 
and defocusing forces from beam space-charge and thermal (omittance) effects. Solutions for rj will be employed to 
seed the IM method of constructing matched envelope solutions. In general, the continuous limit approximations 
tend to be more accurate for weaker applied focusing strengths with aoj ^ 80°. However, even for higher values of 
aoj < 180°, the formulas can still be applied to seed iterative numerical matching methods if the methods have a 
sufficiently large "basin of attraction" to the desired solution. 

For case parameterizations (specified aoj, Q, and £j) the solutions of Eq. I|25|l will, in general, need to be calculated 
numerically from a trial guess. Certain limits are analytically accessible and often relevant. If the beam perveance Q 
is zero, or equivalently if Oj — o"oj, then Eqs. I|25|l decouple and arc trivially solved as 



(ooj/Lp) 



(26) 
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Alternatively, this result can be obtained using Tj = rj in Eq. (|14|) with (Tj — > <7(y . In the case of a symmetric system 



with <j 0x = o- 0y = cr and e x — e y = 
quadratic equation in r? is solved as 



e, then r x 



rf, and the envelope equations decouple and the resulting 



n 



Q , 1 
2 + 2 



41^ 



1/2 



(27) 



In parameterization cases 1-3, the continuous limit solutions TJ must be expressed using the depressed phase 
advances (Xj to eliminate one or more of the parameters Q and Ej. In these cases, if the emittances Ej are known, then 
Eq. lfH|) can be employed to estimate 



((Tj/Lp)' 



(28) 



Alternatively, if the perveance Q is known but one or more of the emittances Ej is unknown, we can use Eq. (|28|l to 



eliminate the emittance term(s) in Eq. 125(1 obtaining (ctq 7 - — (Xj)rj = 2QL^/( 
y-equations yields 



,). Taking the ratio of the x- and 



'Ox 



a 0y ~ °l ' 



(29) 



Back-substitution of this result in (crfy — (Xj)rj = 2QL v j(r x + r y ) then gives 

_ V2QL p 



2] i 

) 



12 1\ I 

(°0x x ) + (a* _oj; 

■\/2QLp 



(30) 



Alternative smooth- limit formulations in the literature |l4l Il5j can also be employed to estimate the rj for systems 
with high degrees of symmetry. 



IV. NUMERICAL ITERATIVE METHOD FOR MATCHED ENVELOPE CALCULATION 

We formulate a numerical iterative matching (IM) method to construct the matched beam envelope functions Tj (s) 
over one lattice period L p using the developments in Sec. II I II The IM method is formulated for arbitrary periodic 
focusing functions Kj. Constraints necessary to apply the IM formalism to all cases of envelope parameterization^ 
listed in Table U are derived. 

Label all quantities varying with iteration number with a superscript i (i = 0, 1, 2, • • • ) denoting the iteration 
order. For example, the ith order envelope functions are labeled r*. The iteration label should not be confused with 
the initial coordinate Si and the initial "seed" iteration corresponds to i = 0. Parameters such as the perveance Q 
or emittances Ej will also be superscripted to cover parameterization cases where the quantities are unspecified and 
are calculated from the envelope functions and other parameters (see Table [|J. For example, e* denotes the j-planc 
emittance at the ith iteration and for parameterization cases where the value of Ej is specified, then = Ej = const. 

For iterations i > 1, we calculate refinements of the principal orbit functions [see Eqs. (0 and I jllj l] in terms of the 
envelope calculated at the previous, i — 1 iteration from 

20 i ~ 1 F i 



j 



Here, denotes Cj(s\si) or Sj(s\si) which are subject to the initial (s = Si) conditions Cj(si|s,) = 1, Cj '(sj|sj) = 
and Sj(si\si) = 0, Sj'(si\si) = 1. Note that the depend on the envelope functions and perveance of the prior, 



8 



i — 1, iteration, previous iteration envelope functions. Updated envelope functions rj and/or betatron functions /3j 
are calculated [see Eq. Q22|l] from the Fj for all i from 



[sj-H*)] 2 



S'j(s l + £ p |si)/ sinaj 
S*( Sl +L p \s 



Sill cr' 



CU Si +L p \ Si ) 



S l Asi + L p \si) 



S]{s\ Si ) 



(32) 



Here, if the parameterization does not specify the depressed phase advances as cr* = Oj, then they are calculated [see 
Eq. O] for all i from 



-[qisi + Lpls^+Sfisi + Lplsi)] 



(33) 



In parameterization cases to 3 (see Tabic HJ), one or more of the needed quantities among Q l , e*, and crj are not 
specified (e.g., for case 1: £) ^ £j specified) and must be calculated to apply Eq. 1(32(1 and/or to calculate the next 
(i + 1) iteration principal functions from Eq. 1)31(1 . Equations 1(33(1 and/or the constraint equations 1(23(1 with Q — > Q l , 
r l j (or in some cases rj — > Js l jP) ) can be employed to calculate parameter eliminations necessary 



, and r,- 



to fully realize each iteration as follows for each case: 

Case (Kj, Q, Ej specified) The cr* can be calculated from Eq. 1(33(1 . 

Case 1 (nj, Q, and Oj specified) The e* can be calculated using Eq. 1(23(1 expressed in betatron form to obtain 



Ik. 

2Q* 



v 

2Q l 



(34) 



with the ratio e y / e l x on the right hand side of the equations determined by 



KxV^-l/(/%) 3/2 



(35) 



Note that expressing the constraints in terms of betatron functions /?j is necessary in this case because the 
envelope functions cannot be calculated from Eq. 1(32(1 until the e* are known, whereas because of the structure 
of the envelope equations, the /3] = (rj) 2 /e* can be calculated from Eq. 1(32(1 without a priori knowledge of the 
values of the . 

Case 2 (Kj, £j, and o~ x specified; or Kj, Ej, and o~ y specified) If necessary, either a l x or o~ y can be calculated 
from Eq. 1(33(1 to enable full specification of the functions /3j or rj. Then, Q % can be calculated using Eq. 1(34(1 
and the f3), or alternatively, using 



- (e'-) 2 /(rj)3 



2Q 4 = 



(36) 



with 



Case 5 (Kj, o~j, and specified; or Kj, <Jj, and e y specified) First, Eq. I(35|) and the functions can be applied 
to calculate e y from specified e x , or e^ from specified e y . Then, Q 1 can be calculated from the e* (if specified, 
= £_,) using Eq. ((34(1 and the ffj, or alternatively, with Eq. ((36(1 and the r*. 
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The seed i = iteration is treated as a special case where the continuous limit formulas derived in Sec. IIIIBI are 
applied to estimate the leading order defocusing effect of space-charge on the beam. In this case the principal functions 
are calculated from 



2QFf 



Ff " + Kj Ff - 7 =^= = 0. (37) 



Here, F® denotes C?(s|sj) or Sj(s\si) subject to the initial (s = Sj) conditions C°(sj|sj) = 1, C?'(sj|sj) = and 
Sj(si\si) = 0, Sj'(si\si) = lj and Q and rj denote the continuous focusing approximation perveance and envelopes 
calculated from the formulation in Sec. IIII 51 with Q — ► Q and Ej — > ej. The continuous focusing values of Q an d 37 
used in calculating the FT are set by the parameterization values in cases where they are specified (e.g., Q — Q for Q 
specified). Otherwise, Q and/or the ~e] arc calculated in terms of other parameters using the appropriate constraint 
equations from Eqs. (|26[I - H30[) applied with Q — ► Q and ej — ► ej. 

Note that the seed envelope functions calculated under this procedure are not the continuous limit functions 
(i.e., r° 7^ rj). Likewise, in parameterizations where they are not held fixed, the seed perveance and emittances will 
not equal the continuous focusing values (i.e., Q° ^ Q and/or £° ^ ~e]). Due to Eq. (|32|l . the seed envelope functions 
r*j will have a (dominant) contribution to the envelope flutter from the applied focusing fields of the lattice with 
a correction due to space-charge defocusing forces from the continuous limit formulas. This approximation should 
produce seed envelope functions r° that are significantly closer to the actual periodic envelope functions rj than would 
be obtained by simply applying continuous limit formulas (i.e., taking r° = rj) or neglecting the effects of space-charge 
altogether [i.e., by calculating r° using Eq. (1221) with Q = 0]. Generally speaking, a seed iteration closer to the desired 
solution can reduce the number of iterations required to achieve tolerance, and more importantly, can help ensure a 
starting point within the basin of attraction of the method, thereby reducing the likelihood of algorithm failure. At 
the expense of greater complexity and less lattice generality, alternative seed iterations can be generated using low 
order terms from analytical perturbation theories for matched envelope solutions 0, 0, S 0- I n certain cases, these 
formulations may generate seed iterations closer to the matched solution. 

Iterations can be terminated at some value of i where the maximum fractional change between the i and (i — 1) 
iterations is less than a specified tolerance tol, i.e., 



Max 



j j 



< tol, (38) 



where Max denotes the maximum taken over the lattice period L p and the component index j = x, y. Many numerical 
methods will be adequate for solving the linear ordinary differential equations for the principal functions Cj and Sj of 
the iteration because they are only required over one lattice period. Generally, the principal functions will be solved 
at (uniformly spaced) discrete points in s over the lattice period. These discretized solutions can then be employed 
with quadrature formulas to calculate any needed integral constraints to affect the envelope parameterizations given 
in Tabled Finally, the convergence criterion l|38fl can be evaluated at the discrete s-values of the numerical solution 
for r l j at the ith iteration using saved i—1 iteration values for r 1 ^ 1 that are needed for calculation of the ith iteration. 
It is usually sufficient to evaluate the convergence criterion at some limited, randomly distributed sample of s-values 
within the lattice period. Issues of convergence rate and the basin of attraction of the method are paramctrically 
analyzed for examples corresponding to typical classes of transport lattices in Sec. 



V. EXAMPLE APPLICATIONS 



In this section we present examples of the IM method developed in Sec. lIVI to construct matched envelope solutions 
and explore parametric convergence properties for example solenoidal and quadrupolc periodic focusing lattices. For 
simplicity, examples are restricted to plane-symmetric focusing lattices with equal undepressed particle phase advances 
in the x- and j/-planes (i.e., o-qj = cto) and a symmetric beam with equal emittances in both planes (ej = e). Under 
these assumptions, the depressed phase advances <7j are also equal in both planes (oj = a) and parametrization cases 
2 and 3 of Table [I] are identical. First, parameterization cases 1 (specified Kj, Q and a) and 2 (specified Kj, e, and 
a) are examined in Sec. IV Al Both of these cases have specified depressed phase advance a and represent the most 
"natural" parameterization of the IM method. Then results in Sec. IV Al are extended in Sec. IV Bl to illustrate how the 
IM method can be applied to parameterization case (specified Kj, Q, and e) with unspecified a while circumventing 
practical implementation difficulties. 
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For simplicity, we further restrict our examples to periodic solcnoidal and quadrupole doublet focusing lattices 
with piecewise constant focusing functions Kj(s) as illustrated in Fig. ^ Solenoidal focusing has k x (s) = K y (s), 
and alternating gradient quadrupole focusing has k x (s) = —k v (s). For both the solenoid and quadrupole lattices 
illustrated, r\ € (0, 1) is the fractional occupancy of the focusing elements in the lattice period L p . The focusing 
strength of the elements is taken to be \nj\ = k = const within the axial extent of the optics and zero outside. 
For solenoids, Kj = k > in the focusing element; and for quadrupoles k x = — K y = k > in the focusing-in-x 
element of the doublet, and k x — ~n y — —k < in the defocusing-in-x element. The free drift between solenoids 
has axial length d = (1 — T))L p . For quadrupole doublet focusing, the two drift distances d\ = a(l — r\)L v and 
&i = (1 — a)(l — r))L p separating focusing and defocusing quadrupoles can be unequal (i.e., d-i). A syncopation 
parameter a G [0, 1] provides a measure of this asymmetry. Without loss of generality, the lattice can always be 
relabeled to take a G [0, 1/2], with a = corresponding to the focusing and defocusing lenses touching each other 
[di = and d.2 = (1 — rf)L p ] and a = 1/2 corresponding to a so-called FODO lattice with equally spaced drifts 
[d\ = d 2 = {l — rj)L p /2\. These focusing lattices are discussed in more detail in Ref. including a description of how 
the focusing strength parameter k is related to magnetic and/or electric fields of physical realizations of the focusing 
elements. 



kJs) 



a) Periodic Solenoid 

( K x = K y ) 




b) Periodic Quadrupole Doublet 
(K x = HC y ) 



FQuad 


dj 
— ► 


riV 2 




d 2 







■nV 2 




DQuad 





Lattice Period 



A 

K 



A 
-K 



dj = a(l-r\)L p 
d 2 = {\-a){\-T\)L p 



FIG. 1: Periodic solenoid and quadrupole focusing lattices with piecewise constant Kj 



As mentioned, for general lattices the scale of the focusing functions Kj can be set by the undepressed phase advances 
&oj using Eq. (|10fl . For the piecewise constant Kj defined in Fig.^ this calculation Q shows that the focusing strength 
|k| is related to CTo by the constraint equations: 

{cos(26) - sin(29), Solcnoidal Focusing, 

cos 6 cosh 6 Quadrupole Focusing. 

+ ^20(cosesinhe -sinOcoshe) ( 39 ) 
- 2a(l-a)ii^9 2 sinesinhe. 

Here, for both solenoidal and quadrupole focusing lattices, 9 = y/\k\rjL p /2. In the analysis that follows, Eq. (|39ll is 
employed to numerically calculate for a specified value of <7o, and then k is calculated in terms of other specified 
lattice parameters as \k\ = 49 2 / (r)L p ) 2 . The undepressed phase advance cto is measured in degrees per lattice period. 
Integrations of needed principal orbits to implement the IM method are carried out with the with initial conditions 
(s = Si) corresponding to the axial middle of drifts separating focusing elements. 

Typical matched envelope solutions tj are shown for one lattice period in Fig. [3 for solenoid, FODO (a = 1/2) 
quadrupole, and syncopated (a 7^ 1/2) quadrupole focusing lattices. Scaled x-plane lattice focusing functions k x are 
shown superimposed. Excursions of the matched envelope functions are in-phase for solcnoidal focusing (r x = r y ) 
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because the applied focusing is plane symmetric (k x = n y ). In contrast, for quadrupolc focusing, the anti-symmetric 
plane focusing (k x = —K y ) results in out of phase envelope flutter in each plane (focus-defocus) leading to net 
focusing over the lattice period in both planes. Expected symmetries of the matched solutions are present for both 
the solenoidal and quadrupole focusing lattices (see Appendix |A"|) . For the quadrupole solutions, note that the FODO 
case exhibits a higher degree of sub-period symmetry than the syncopated case. Leading order terms of an analytical 
perturbation theory for the matched beam envelope solutionQ can be applied in the limit a — > to show that the 
envelope excursions (flutter) scales as 



MaxL 



(l-cos g )(l-??)(l-'>7/2) 
6 

(l-cos g ) 1 / 2 (l-r ) /2) 
2 3 /2(i_2) ) /3) 1 /2 



Solenoidal Focusing. 
FODO Quadrupole Focusing 



(40) 



Equation (|4U|) shows that for solenoidal focusing the matched envelope flutter increases with decreasing lattice occu- 
pancy r\ and increasing focusing strength ero. In contrast, for FODO quadrupole focusing the flutter depends weakly 
on r] (the variation of Max[r x ]/r^ — 1 in rj has a maximum range of 0.07) and more strongly on a (variation of 0.5). 
Envelope flutter changes only weakly when space-charge strength is reduced (i.e., cr/a increased). 




0.4 0.6 0. 

s/L p 



FIG. 2: (Color) Example matched envelope solutions for (a) solenoidal, (b) FODO (a — 1/2) quadrupole, and (c) syncopated 
(a = 0.1) quadrupole focusing lattice. Parameters for all cases are: L p = 0.5 m, 77 = 0.5, an = 80°, Q = 4 x 10" 4 , and e = 50 
mm-mrad. These parameters yield a/ao — 0.3144, 0.3093, and 0.3099 for (a), (b), and (c). 



Although the system symmetries assumed simplify interpretation of the matched envelope solutions obtained in the 
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examples, we note that the numerical methods employed in calculation of the particle principal orbits functions and 
any necessary constraint equations are not structured to take advantage of the symmetries of the matched solutions. 
Because of this, the examples provide a better guide to the performance of the IM method in situations where there are 
lesser degrees of system symmetry. Mathematical^ based programs used to in the examples have been archived [l7j. 
These programs can be easily adapted to more complicated lattices. 



A. Case 1 and 2 parameterizations 

The IM method described in Sec. IIVI is applied with a 1 = a specified and the unknown parameters of the ith 
iteration e % (case 1: Q and a specified) or Q l (case 2: e and a specified) calculated from the constraint equations 134(1 
and 1(35(1 . The continuous focusing approximation envelope radii TJ used in the seed (i = 0) iterations are calculated 
from Eq. I(3U|I (case 1) and Eq. 1 (28(1 (case 2). The number of iterations needed to achieve a 10 -6 fractional envelope 
tolerance [see Eq. ffify ] are presented in Fig. as a function of cr an d cr/cro for solenoidal, FODO quadrupole, and 
syncopated quadrupole focusing lattices employing both case 1 and case 2 parameterization methods. In Fig. 0] 
iterations corresponding to one data point in Fig. [3] (case 2 solenoidal focusing) are shown. This example shows a 
result typical for lower to intermediate values of <tq , where the seed iteration r° is fairly close to the matched solution 
and the first iteration correction rj closely tracks the matched solution to within a percent local fractional error. 
Higher values of and more complicated lattices result in both seed iterations farther from the matched solution 
and less rapid convergence with iteration number. 
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c) Case 2, e 
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FIG. 3: Number of IM iterations needed to achieve a tol = 10~ 6 fractional error tolerance matched envelope solution for (a) 
solenoidal, (b) FODO (a = 1/2) quadrupole, and (c) syncopated (a = 0.1) quadrupole focusing lattices as a function of ao (for 
(Jo = 40°, 60°, 80°, ■ ■ ■ , 160°) and a/ao (for a/ao = 0.1, 0.2, 0.3, • ■ ■ , 1.0). The left column corresponds to the parametrization 
case 1 method with Q = 10~ 4 (Unachievable limit points marked x) and the right column corresponds to the parameterization 
case 2 method with e = 50 mm-mrad. Other lattice parameters are L p = 0.5 m and 77 = 0.5 for (a), (b), and (c). 
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0.0 0.2 0.4 0.6 0.8 1.0 

s/L p 

FIG. 4: (Color) For solenoidal focusing: converged matched envelope solution (black) with r x = r y , the first seed iteration 
(red), and the second iteration (green). System parameters are L p = 0.5 m, 77 = 0.5, cto = 80°, e = 50 mm-mrad, and 
a /<jo = 0.3. 

The data in Fig. |3| shows that the IM method converges rapidly to small tolerances over a broad range of ap- 
plied focusing (<7o) and space-charge (ct/cto) strength. Not surprisingly, stronger focusing strength (i.e., increasing 
cto) requires more iterations for both solenoidal and quadrupole focusing at the fixed value of lattice occupancy 77 
employed. Also, lesser degrees of lattice symmetry result in more iterations being necessary for convergence (e.g., 
lattice convergence rate order: solenoidal, FODO quadrupole, syncopated quadrupole). Iterations required appear 
to depend only weakly on space-charge strength (ct/ct ) - except for solenoidal focusing lattices with very high ct 
where required iterations become abruptly larger for weak space-charge with ct/cto close to unity. Even parameters 
deep within the regime of strong envelope instability Q converge rapidly. Points for ct/cto = 1 arc eliminated in the 
case 1 examples because the perveance Q is held to a fixed, finite value and this limit would correspond to a matched 
beam envelope with infinite cross-sectional area. Conversely, for the limit ct/ctq = 1 in the case 2 examples, only one 
iteration is required for convergence to a finite solution because for zero space-charge strength the trial seed iteration 
generated by Eq. l).'i2l) for i = corresponds to the exact matched envelope to numerical error (i.e., when Q Q = 
the C® and S® are the principal undepressed particle orbits which generate the matched envelope of the undepressed 
beam). The IM method applies for extremely strong space charge with ct/cto <C 0.1, but probing the limit a — > 
requires careful analysis of various terms in the formulation presented in Sec. IIVI 

Complementary to Fig. [3J the decrease in the log of the fractional tolerance [see Eq. 1|38[) ] achieved with iteration 
number is plotted in Fig. [S]for solenoidal and FODO quadrupole focusing lattices for one set of system parameters. 
The matched envelopes are calculated using the case 2 methods. Case 1 methods and other system parameters 
yield similar results to those presented. We find that the IM method converges rapidly, with the fractional tolerance 
achieved increasing by one to two orders of magnitude per iteration till saturating at a value reflecting the precision 
of numerical calculations employed (~ 10~ 15 fractional accuracy for the examples in Fig. [SJ. 
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FIG. 5: Log of fractional error tolerance (tol) achieved for matched envelope solutions versus number of IM iterations for (a) 
solenoidal and (b) FODO (a = 1/2) quadrupole focusing lattices. Solutions are generated using case 2 methods for system 
parameters: L p = 0.5 m, r\ — 0.5, tro = 80°, e — 50 mm-mrad, and a/ao = 0.2 [corresponding to Q ~ 6.700 x 10~ 4 and 
Q ~ 6.561 x 10~ 4 for (a) and (b)]. 



B. Case parameterization 

In parameterization case 0, the matched envelope functions rj are specified by Kj, Q, and Ej. For the ith iteration, 
the depressed phase advances <r*- needed to calculate the iteration envelope functions r* are most simply calculated 
using Eq. Unfortunately, this simple method can fail if space-charge is strong and iterations result in envelope 

corrections where the radial cross-section of the beam is compressed sufficiently relative to the actual matched envelope 
solution over the lattice period. Compressive over-corrections can produce iteration principal orbits that are depressed 
below zero phase advance (i.e., er* =0). In this situation, er* becomes complex and we find that Eq. (|32f> can fail to 
generate an iteration closer to the desired matched envelope solution. 

To better understand where the problem described above can occur, a simple continuous focusing estimate (see 
Sec. IIII B"|) is applied. Taking Kj = (ao/L p ) 2 and Ej = e, we estimate the envelope compression factor / needed to 
fully depress particle orbits within the matched envelope. A particle moving within the continuous matched envelope 
Tj = Tb has depressed phase advance (a/L p ) 2 = (ao/Lp) 2 — Q/1% 2 . Replacing j% — ► ff% and a — > in this phase 
advance formula gives 

* (41) 



But for continuous focusing, we havejj| 



^1+ ^T+Ma e/(QL p )] 



aoe_ (0-/00) ( , 9] 
QL P I-Wct-o) 2 " K ' 



Together, Eqs. and show that / = 0.99, 0.95, and 0.90 (corresponding to ~ 1%, 5% and 10% compressive 
over-corrections) will produce fully depressed particle orbits for cr/cro < 0.14, 0.31, and 0.44. Numerically analyzed 
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examples below indicate that this problem can occur in periodic focusing lattices for more moderate space-charge and 
compression factors than the continuous focusing estimates suggest. 

The parameter region where the IM method can be applied using the "conventional" case procedure for example 
periodic solenoid and FODO quadrupole lattices is illustrated in Fig. |SJ The region of applicability corresponds to 
parameters where Eq. (|33[l can be employed to calculate the iteration depressed phases advances cr* without obtaining 
complex values. Iterations necessary to achieve tolerance are plotted as a function of ao and <j/<7q. Rather than 
plotting results in terms of the perveance Q, Eq. I|14f> was used to calculate <j/<jq from the matched envelope functions 
and system parameters to better quantify the relative space-charge strength where the method fails. Values of Q were 
chosen to uniformly distribute points in cr/uo- Note that the IM method works with the simple initial seed iteration 
when space-charge is moderate to weak (0.6 < a/ao < 1) but abruptly fails with increasing space charge (o-/tr < 0.6). 
Near the point of failure, convergence becomes slow (iteration counts for the example in Fig. [5] can become thousands 
if points are chosen sufficiently close to the start of the failure region) . 
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FIG. 6: Number of "conventional" IM iterations needed in case to achieve a tol = 1CF 6 fractional error tolerance matched 
envelope solution for (a) solenoidal and (b) FODO (a — 1/2) quadrupole lattices as a function of ao (for ao = 40°, 60°, 80°, 
• • ■ , 160°) and a/ao (for a/ao = 0.1, 0.2, 0.3, • ■ ■ , 1.0). System parameters are: L p — 0.5 m, 77 = 0.5, and e — 50 mm-mrad. 
Parameters where the method fails are marked x. 

Several alternative methods were attempted to render the IM method applicable to all case parameters with 
arbitrary space-charge strength. We describe these methods without assuming aox = &0y and e x = e y to better reflect 
general case applications. First, rather than employing Eq. to calculate the depressed phase advance <rj of the 
iteration, the integral formula (|14f) is applied with the envelope functions of the previous i — 1 iteration with 

<r;=e,-/ . (43) 

The anticipation is that cr* calculated from Eq. I|43() should be sufficiently close to the actual depressed phase advance 
Uj of the converged solution to correct the problem. Unfortunately, this method, when applied to example solenoid 
and quadrupole lattices, results in systematic convergence to unphysical solutions. Replacing Eq. I|43|) with an "under- 
relaxed" average over previous iterations might address this problem but was not analyzed. In cases where complex 
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phase advances resulted, various other simple replacements of Eq. I|33|) were attempted without obtaining satisfactory 
results. 

Several alternative procedures extend applicability to general case parameters. First, slowly increasing the per- 
veance Q from some sufficiently small (or zero) value while implementing the conventional case iteration method 
using Eq. I|33[) proves workable in our tests. In this scheme, if Eq. (|33() fails (i.e., produces unphysical complex values 
for cjJ) then Q is adaptively decreased while iterating until the formula becomes valid before increasing Q again 
toward the target value. For strong space-charge this procedure can result in many iterations being necessary for 
convergence because small increases in Q were required in various test cases examined. It is also difficult to determine 
optimal increments to increase the pcrvcance - which complicates practical code development and can limit the range 
of method applicability. 

Another, simpler to implement, alternative procedure is formulated by combining the Sec. IV Al method for solving 
case 1 parameterizations with numerical root finding. In this "hybrid" procedure, the emittances Ej calculated from 
the x- and y-plane constraint equations 1|23[1 are regarded as an undetermined function of the Oj [i.e., £j| S pecified = 
£j [cr x , (J y )] and trial matched envelope solutions Vj are rapidly calculated to tolerance using matched envelopes obtained 
with case 1 methods for specified (guessed) values of the aj. Numerical root finding can be employed to refine the 
guessed values for the <jj to obtain the values of Oj consistent with the target values of Ej. Because the Ej(a x , a y ) are 
smooth, monotonic functions of the Oj for < cry < aoj > the consistent values of the <jj can be found with relatively 
small numbers of root finding iterations. This is particularly true for plane-symmetric systems (troj = Co and Ej = s) 
because one-dimensional root finding can be employed. 

The total number of two-dimensional (i.e, the calculations do not assume plane symmetry) iterations needed to 
implement this hybrid method for case is shown in Fig. [7| for example periodic solenoid and FODO quadrupolc 
lattices. Here, the total iteration number represents the sum of all iterations needed to calculate the emittances to 
a specified fractional tolerance while calculating all trial matched envelope solutions to a separate specified tolerance 
over all two-dimensional root finding steps. The same lattices and presentation methods used in Fig. are employed 
to aid comparisons. Note that the full case parameter space is accessible in this procedure with only relatively 
modest total iteration counts in spite of the additional numerical work resulting from the root finding. A secant-like 
multi-dimensional root finding method is emploved|l6j. Note that only two-dimensional root finding is necessary in 
contrast to four-dimensional root finding associated with conventional procedures for constructing matched envelope 
solutions by finding appropriate initial envelope coordinates and angles. Initial root finding iterations are seeded using 
continuous focusing model estimates for aj calculated from Eq. (|28[1 using the seed values of rj. Subsequent root 
finding steps in <jj employ the previous step matched envelopes as a seed envelope in the case 1 iterations. For small 
root finding steps in Oj this previous step seeding saves considerable numerical work. Only one iteration is necessary 
for the limit points with ct/ctq = 1 because for zero space-charge strength the trial seed iteration is exact to numerical 
error. Iteration counts at fixed ctq likely increase and decrease in <r/cr due to approximate iteration seed guesses being 
(accidentally) farther and closer to the actual root than in other cases. If the plane symmetries are employed (i.e., 
using e x = £ y and a x = <j y ), then total iterations required can be further reduced. Matched envelopes for general case 
parameters can also be calculated in similar number of total iterations by analogously combining case 2 methods 
with numerical root finding. In this case values of <jj consistent with specified values of Q are calculated using the 
two components of Eq. 112,'jH [i.e., set Q — > Qj in the j = x and y components of Eq. 112.'j[) and then root solve for Gj 
consistent with Q = Qj(a x ,a y )]. 
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FIG. 7: Number of total "hybrid" IM iterations needed for case using tol = 10 fractional error tolerance case 1 matched 
envelope solutions and root found £j with 1CP 4 fractional accuracy from specified values for (a) solenoidal and (b) FODO 
(a = 1/2) quadrupole lattices. The same system parameters and presentation format is employed as is used in Fig. [fj] 



VI. CONCLUSIONS 

An iterative matching (IM) method for numerical calculation of the matched beam envelope solutions to the KV 
equations has been developed. The method is based on orbit consistency conditions between depressed particle 
orbits within a KV beam distribution and the envelope of orbits making up the distribution. Application of the 
IM method in simplest form requires numerical solution of linear ordinary differential equations describing principal 
particle orbits over one lattice period and the calculation of a few axillary integrals over the lattice period. A large 
basin of convergence enables seeding of the iterations with a simple trial solution that takes into account both the 
envelope flutter driven by the applied focusing lattice and leading-order space-charge defocusing forces. All cases of 
envelope parameterizations can be employed, but the method is most naturally expressed, and highly convergent, when 
employing the depressed particle phase advances <jj as parameters — which also corresponds to a natural choice of 
parameters to employ for enhanced physics understanding. Virtues of the IM method are: it is straightforward to code 
and applicable to periodic focusing lattices of arbitrary complexity; it is efficient for arbitrary space-charge intensity; 
and it works for all physically achievable system parameters - even in bands of parametric envelope instability where 
conventional matching procedures can fail. 
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APPENDIX A: MATCHED ENVELOPE SYMMETRIES FOR QUADRUPOLE DOUBLET AND 

SOLENOID AL FOCUSING 

Consider a periodic quadrupole doublet lattice 0] focusing a beam with symmetric emittances (i.e., e x = e y ). To 
concretely define doublet focusing, we assume that the lattice focusing functions Kj(s) satisfy 

Kj(s) = -Kj(-s), (Al) 

in addition to the general quadrupole lattice symmetry n x = —K y . This doublet focusing symmetry is consistent with 
focusing/defocusing elements with axial structure (i.e., including fringe fields) if both the focusing and defocusing 
elements are realized by identical hardware assemblies with equal field excitations appropriately arranged in a regular 
lattice via symmetry operations (i.e., translations and rotations). Without loss of generality, let s = correspond 
to the axial location of the drift between two successive quadrupoles in the periodic lattice (for cases where a finite 
fringe field extends into the drifts, this location will be where Kj = 0). Assume that the matched envelope functions 
satisfying the KV equations Q arc symmetric about the mid-drift with 

r j (s)=r- j (-s). (A2) 

Here, if j = x,y, then j = y,x. Take the j = x KV equation [see Eq. Q], substitute s — ► — s. Then employing 
the focusing and envelope symmetries in Eqs. (|A1|) and 1A2|1 together with n y = —k x obtains the complementary 
j = y KV equation, thereby showing that the assumed symmetry in Eq. 1|A2|) is consistent. An immediate corollary of 
Eq. i|A2fl is that at any mid-drift between quadrupoles, the envelope is round (i.e., r x = r y ) with opposite convergence 
angles (i.e., r' x = —r' y ). 

Restrict the situation described above to a symmetric FODO system where the two focusing and defocusing 
quadrupoles are separated by equal length axial drifts^ and the focusing and defocusing elements are each re- 
flection symmetric about their axial midplane [i.e., within one element, k x (s — s) = k x (—s + s) where s = s is the 
geometric field center of the element] . These further assumptions lead to the additional FODO focusing symmetry 

K j (s)=K~ j (L p /2 + s). (A3) 

With the choice of s = made as above, the focusing and defocusing optical elements are centered at s = L p /A and 
s = 3L p /A within the period s £ [0, L p }. Using steps analogous to those outlined above, it can be shown that the 
matched envelope functions also have the FODO symmetry: 

r j (s)=r~ j (L p /2 + s). (A4) 

Another FODO symmetry can be obtained by replacing s — > —s in Eq. l|A4f) . applying Eq. i|A2fl . and differentiating 
to yield 

r' j (s) = -r' j {L p /2-s). (A5) 

Evaluating this expression at the focusing element centers at s = L p /A and s = 3L p /4 and invoking periodicity of the 
Tj with r'j(s + L p ) = r^-(s) shows that the matched envelope functions are extremized (i.e., r'j = 0) at the focusing 
element centers in a symmetric FODO lattice. The envelope equations then shows that the j-plane extrema of 
Tj with Kj < (defocusing plane) satisfies r'j > and therefore must be a minimum value. Period symmetries then 
require that the other focusing plane extrema (j-plane with Kj > 0) corresponds to a maximum value. 

Analogous steps to those employed in the analysis of quadrupole doublet focusing can be applied to solenoidal 
focusing (k x = K y ) systems with e x — e y to show that r x = r y . Consider a periodic solenoidal focusing function 
with only a single clement in the period that is also reflection symmetric about the axial midplane (with reflection 
symmetry defined as for the FODO quadrupole case above). Then procedures used above are readily employed to 
show that the matched envelope function Tj is maximum at the axial center of the focusing element and is minimum 
at the axial center of the drift. 
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